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Abstract 

The effective formulas reducing the two-dimensional Hermite polynomials to the 
classical (one-dimensional) orthogonal polynomials are given. New one-parameter gen- 
erating functions for these polynomials are derived. Asymptotical formulas for large 
values of indices are found. The applications to the squeezed one-mode states and to 
the time-dependent quantum harmonic oscillator are considered. 



1 Introduction 



The Hermite polynomials of several variables arise quite naturally almost in all problems 
relating to quantum systems described by means of multidimensional quadratic Hamiltoni- 
ans: see, e.g., [1-4] and references therein. However, until now they were not widely used 
by other authors in the papers on quantum mechanics and quantum optics (for a few of 
exceptions see, e.g., [5,6]), due to the absence of simple explicit formulas, which would be 
convenient for the calculations and for the analysis of the relations obtained. The aim of 
the present paper is to give the expressions for multidimensional Hermite polynomials in 
terms of the well known classical orthogonal polynomials. Besides, we shall give new sum 
rules and generating functions, as well as asymptotic formulas for various combinations of 
the parameters. Although some of the results could be found in odd form in other references 
[1-11], we beheve that bringing them together will be useful for the further applications. 
Just having in mind these applications, we pay a special attention to the specific sets of pa- 
rameters defining the polynomials, which are typical for the quantum mechanical problems. 
We shall consider mainly the case of the Hermite polynomials of two variables, but some 
generalizations to higher dimensions will be also given. 

The structure of the paper is as follows. In the next section we derive the explicit ex- 
pressions for the two-dimensional Hermite polynomials in terms of the Jacobi, Gegenbauer, 
Legendre, and usual Hermite polynomials. The asymptotics of the two-dimensional Hermite 
polynomials of zero arguments are considered in Section 3. New sum rules and generating 
functions for the "diagonal" two-dimensional Hermite polynomials are given in Section 4, and 
their generalizations to higher dimensions are considered in Section 5. Section 6 is devoted to 
the physical applications. In that section we discuss briefly two problems: the photon statis- 
tics in squeezed mixed (thermal) quantum states, and the transition probabilities between 
the energy levels of a quantum oscillator with a time-dependent frequency. 

2 Relations between two-dimensional Hermite polyno- 
mials and classical orthogonal polynomials of one 
variable 

The two-dimensional Hermite polynomials i7^^^(7/i, 7/2) are defined by means of the gener- 
ating function [12], 
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^ n,m=0 

Here ai and 02 are arbitrary complex numbers combined into the two-dimensional vector 
a= (ai,a2) : 

2 2 
aRa = ttiRikak, aRy = ^ aiRikVk, 

i,k=0 i,k=0 

and R is the symmetric matrix 
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Let us begin with the case when y = (in the quantum mechanical apphcations it 
corresponds, for example, to the absence of an external force). Suppose that Ru ^ and 
R22 7^ 0. Designating the left-hand side of eq. (1) as 



E = 



exp (^—Riial - ^R22al - ^120102 ) , 



we expand it in the following Taylor's series. 
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Introducing the notation 
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we can rewrite the square bracket in the right-hand side 
1 2ra + o? 



as 



\a + rf - (r' - 1)] = (1 - r') [l - (p + 7) 



where 



7 = 



Developing the function [1 — (p -|- 7)^]^ with respect to variable 7 we arrive at the expansion 
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Taking into account the definition of the associated Legendre functions [12], 



2iq\ ^ ' 



we obtain the following formula for the two-dimensional Hermite polynomial of 
ments: 
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(2) 
zero argu- 
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where 



fJ'mn = min(m,n), 

and integers m, n must have the same parity (otherwise the right-hand side equals zero). 
We draw an attention to the definition of the associated Legendre functions: eq. (2) means 



that we consider them in the complex domain, which does not include the segment [—1, 1]. 
Consequently, the definition adopted in the present paper differs from the definition used, 
e.g., in the theory of spherical harmonics and in Refs. [2,3,7], where a special case of eq. 
(3) was derived in connection with the calcTilations of the transition probabilities between 
the energy levels of a quantum oscillator with a timc-dcpcndcnt frequency. The reason is 
that in the generic case we can hardly expect the variable p to belong to the segment [—1, 1] 
(although in some specific but important cases this is possible: see the last section of the 
paper) . 

For coinciding indices we get 

i/W(0, 0) = n!(-detR)tp„ (^^g=) , (4) 

Pn{z) being the usual Legendre polynomial. Using the relations between the associated Leg- 
endre functions P^{z), the Jacobi polynomials P^'^\z)^ and the Gegenbauer polynomials 
C.iz) [12,13], 

one can write eq. (3) in the following equivalent forms. 
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These expressions hold for arbitrary values of the variable p. 

For nonzero vector y function Hj^^{yi, 1/2) can be written as a finite sum of products of 
the usual Hermite polynomials [8], 

.'f^f 2i?i2 V n\m\ ( \h ( \ (7) 

where 



Ci = RiiUi + Ri2y2, C2 = Ri2yi + R22y2- (8) 

The proof is straightforward. First we express the left-hand side of eq. (1) as the product 
of three exponential functions, 

2i?i2 \ 



exp ^— -aRa + aRy^ = exp (—^1 + "^Zirjij exp (^—zl + 222^2) exp ^- 
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where 



/ Ru / R22 Ci C2 

Then we expand each of three exponentials in the power series of Zi, Z2, and Z1Z2, respectively, 
having in mind that the first two exponentials are the generating functions for the usual 
Hermite polynomials. Combining the terms with the powers z^"'^ and -22"^^ from the first 
two expansions with the terms {Z1Z2Y arising in the expansion of the third exponential, 
/c = 0, 1, . . . , n, one arrives at formula (7). Note that it remains valid even when Ru = or 
R22 — 0, provided one takes the corresponding limit. For instance, if R22 — 0, then 



Hi^Hvu y.) = R,2 E [-y) i^n-mm-m\ ^-' [ J " 

Proceeding to the limit Ru = in eq. (9), we get the formula for the two-dimensional 
Hermite polynomial with matrix R = ta^i where Ux is the standard Pauli matrix, 

^-=(1 o)- 

It reads 

Taking into account the definition of the associated Laguerre polynomials [12], 
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we can rewrite (10) as follows. 
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Ht^Hyi,y2) = fimn\t''-{-ir-y, ' y2 ' Ll-7l(iz/iZ/2), (11) 

where 

Urnn = max(m,n). 

For coinciding indices we get 

Hi':-Hy„y2) = (-l)"rn!L„(tyit/2). (12) 

The special cases of eqs. (11) and (12) were obtained earlier in [9,11]. 
Putting y = in eq. (7) we get a formula for i7^^^(0, 0) which may be the most suitable 
for the numerical calculations: 

.2,(0,0) ^ (fip)' t (13) 



where m and n must have the same parity. For m = n we have 



Immediately from the generating function (1) wc obtain the expression for Hj^{yi,y2) 
in terms of iJ|^^(0, 0) and variables Ci, C2 defined by cq. (8): 

Hj^lHyu y2) = E E " 7 ^^^^^^^ o)cr'cr'- (is) 
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3 Asymptotic formulas for i/W(o,0) 

The asymptotics of the functions Hj^^^O,0) for large values of indices can be derived, due 
to eq. (6), from the asymptotics of the Gegenbauer polynomials. Consider first the case 
when the argument of these polynomials, r/\/r'^ — 1, does not belong to the interval [—1, 1] 
(remind that r = -R12/V-R1 1-^22 )■ Then the generalized Laplace- Heine formula (see Theorem 
8.21.10 from ref. [6]) yields 



where 



z^l{C + C'), ICI>1- (17) 
The solution of eq. (17) with respect to ( reads 



^^^ + ev/^^=^^, e = ±l, (18) 
V?^ — 1 

where the choice of sign of e is determined by the requirement |C| > 1. Putting expression 
(18) into eqs. (6) and (16) we get for n » 1, |n — m| ~ 0{1), 

I -t \ "+m \n-m.\ + l 

7/iS'(0,0) « 'J- ' r (!i±f±i) «^ (r + e)'-^. (19) 

We have used the identity [12] 

r{z)T{z + ^) = V^T{2z)2'-'\ 

If parameters R12 and \/Ri\R2i are real, then e = sign(r), and eq. (19) can be rewritten as 
follows, 

n-j-m-|-l 

where 



_ r if < 



if i?12>0 



However, eqs. (19) and (20) are not valid for small values of |r|, since for R12 — and 
odd values of n and m they do not lead to the obvious relation H}^{0, 0) = 0, which is an 
immediate consequence of definition (1). In this case we may use the asymptotic formula for 
the associated Legendre functions given in ref. [14], 



Vshe)-{n+-j 1-^1'/. 
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a, (21) 



n>l, m~C)(l), Re^>0, |Im^| < tt, 

Im{z) being the modified Bessel function (Bessel function of the third kind). Due to eq. 
(18), 

i = logc = log 

Taking into account the relation 



pm/ N _ (n + m)! 
and a consequence of Stirling's formula, 

r(x + y) fa r(x)x^, x>i, y~e)(i), (22) 

we arrive at the expression 



where 
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(23) 



For real r > 1 eq. (23) turns into (20) provided one replaces the modified Bessel function 
by its asymptotic expression for a large argument [15], 



V27fi' 

Now suppose that r is real, but \r\ < 1. Then 



;2| > 1, a -0(1). 



C = ^log ^^^^j - i|e, e = sign(r), 

and we may rewrite eq. (23) in terms of the usual Bessel functions, according to the relation 

Ik{z) ^i-^Jk{iz). 

It is sufficient to make calculations for r > 0, having in mind that the change of sign of r 
results in an extra factor (—1)" = (— l)"* before H^{{), 0), due to eq. (3). Thus we get 
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Taking into account the asymptotics of the Bessel function [12,15], 

Jk{z) ~ — cos — -kn — -71^ , |arg(2;)| < tt, \z\ 3> 1, 
we can simplify eq. (24) as follows. 
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l+|r|)-^'"" + (-l)'*(l-|r|)''"'" . (25) 



If r = 0, this expression exactly equals zero for odd values of n and m (remind that m and 
n must have the same parity). Moreover, it coincides with eq. (20) provided (1 + |r|)'^"'" » 
(1- |r|)-^"'". 

The exceptional case when coefficient p (the argument of the Gegcnbauer or Jacobi 
polynomials) belongs to the interval [—1,1], so that p = cos6, corresponds to r = icot^^. 
Then we use the Hilb-like asymptotics of the Jacobi polynomial, given by Theorem 8.21.12 
ofRef. [13], 



where N — n + {a + P + l)/2, a > —1, and P may be arbitrary number. Putting this 
expression into eq. (5) we arrive at the formula 



H^^HO, 0) - \ Jk^ , (26) 

which is vahd for ^ 1, |ri — m| ~ < < vr. 

Eq. (6) enables also to find the asyptotics of Hj^^^O, 0), when \m — n\ 1, Unm C'(l). 
The definition of the Gegenbauer polynomial [12], 

with account of eq. (22) , yields for A » 1 

[f] /'_1^A;\g-fc 

Comparing this expression with the definition of the Hermite polynomial [12], 

we get the asymptotic formula for the Gegenbauer polynomial, 

Ciiz)^^H,{z^/X), A»l, g~0(l), (27) 
and its consequence for the two-dimensional Hermite polynomial of zero arguments. 
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n-m|>l, //„^~(!?(1). (28) 



(We have used again eq. (22)). The hmit cases of formula (28) for small and large values of 
the argument of the Hermite polynomial are quite clear. 



4 One-parameter generating function and sum rules 
for the "diagonal" polynomials 

Generating function (1) has two auxiliary parameters, ai and 02. However, for the "diagonal" 
polynomials, with m — n, a one-parameter generating function exists. To find it we choose 
vector a as follows. 



where A is a new complex parameter. Then multiplying both sides of eq. (1) by the factor 
exp(— and calculating the gaussian integral 



n,m=0 

with the account of the relation 
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G{X) = / exp(- JaRa + af - (3(3*)d(3d(3* 

J p*''p"'e-^f^*dpdp*, 
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we arrive at the formula (I means the unity matrix) 
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= E HR-"'f) = [det (Aa.R + I)]"^ exp 



The explicit form of the right-hand side of eq. (29) reads (f = (/i, /2)) 



(29) 



G(A) 



1 + 2Ai?i2 - A^ det R 



■I / 2A/1/2 - A^ (/2i?22 + fiRu - 2M2R12) 

exp < 



2 [1 + 2Ai?i2 - A2 det R] 



(30) 



If f = 0, then this function is nothing but the generating function of the Legendre polynomials 
[12], and we again arrive at eq. (4). 

Eq. (30) enables us to represent the "diagonal" two-dimensional Hermite polynomials in 
the form of a finite sum over the products of the Laguerre polynomials. For this purpose we 
rewrite the right-hand side of eq. (30) as follows, 
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exp 
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where 



\i — {'\[RiiR22 — R12) ^ I A2 — — (y^ R11R22 + R12) ^, 

flR22 + - 2/F?^/i/2 
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4-\/ R\\R22 (V R11R22 — -R12) 

/1-R22 + /2-R11 + 2V-R11-R22/1/2 



4a/ R11R22 {\/ R11R22 + R12 
Using the known relation [12] for the Laguerre polynomials 



we see that the generating function G{X) can be represented as 
Thus we obtain the formula 



i/W(R-if) = (-l)"n! J2{R,2 - V^iii?22)'(i?i2 + ^RnR22r-'Ls ^x,)Lj,{x2). (31) 

Assuming f = Ry, we get quite the same decomposition (31) for the function H^^^{y), but 
with the variables Xi and X2 expressed directly in terms of yi and 1/2 as follows, 

2^1 = ^ ^1 - ^^l^^^j + 1/2^22 - 2^i?iii?22|/l2/2) , 

X2 = i (^1 + (2/1^11 + 1/2^22 + 2^RnR22yiy2^ . 

Formulas like (30) and (31) for the diagonal matrix elements of the Gaussian density 
matrix in the Fock basis were obtained also in Rcf. [16]. 

Differentiating the generating function G{X) k times with respect to A and putting A = 1 
in the final expression, one can calculate the cumulants < n{n — 1) • • • (n — /c + 1) >, defined 
according to the relation 

£ n(n - 1) • • • (n - A; + i)l^m!M. ^ ^ ri{n - 1) ■ ■ ■ {n - k + 1) > H^nHy) _ 

n=0 ^- n=0 ^■ 

The first two derivatives yield the "average value" < n > and the "dispersion" 

(7„ =< > - < n >^=< n{n -l)> + <n>-<n>^. 
Differentiating the preexponential factor in eq. (29) with the aid of the matrix formula 

|logdetAW^T.(A-f), 

we get the following general concise relations, 

< n > (R, y) = ^TrS - 1 + ^yRSV.Ry, (32) 

a„(R, y) = iTr[S(S - I)] + iyR(2S - I)SV,Ry, (33) 
where matrix S is defined as follows, 

S = ((7,R + I)-\ 

Of course, these formulas hold provided the series converge. The simplest special case 
corresponds to the matrix R = ta^, when 

t f t V t f^C^-t) 

<'^>= -^r—L + y^y^i^r—:] > = -71— a? + 1/11/2- 



l + t ' "'"'\l + tj ' " (l + i)2 ' ^^^^(i + i)3- 
These relations have sense for — 1 < i < 0. 



5 Generalizations to higher dimensions 

To generalize the results of the previous section to the "diagonal" Hcrmite polynomials of 2N 
variables we introduce complex parameters Aj, i — 1,2, . . . , N, and choose the auxiliary 
parameters aj as follows, 



Then the same scheme results in the following expansion, 
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= [det (AE^ + l2jv)] ^exp 



JyR(AE,R + l2iv)~'S,ARy 
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(34) 



Here y = {yi,y2, ■■■y2N), the 2Nx2N matrix E^; is the 2A'"-dimensional analog of the Pauli 
matrix ax. 



" ' Im ' 



In is N X N unit matrix, and the diagonal 2Nx.2N- matrix A reads 

N 

where each matrix Aj has only two nonzero elements: 

= i^j)j+N,j+N = 1- 

Note that matrix A commutes with matrix E^. but does not commute, in general, with matrix 
R. 

Generating function (34) enables to calculate the "average values" < ni^ri^^ ■ ■ -n^^ >, 
defined as follows. 



OO OO 



J2 ■■■ 12 W^'n^-'-n}, 



if{R-} (v) 
ni=o n^=o ni!n2!---niv! 

OO OO H-{R.} fvl 

<- III 112 ' " ''■N ^ / . ■ ■ ■ / . 



nT^o nT^o ni\n2\---nN\ 
Taking into account the relations 
dA 

— = Aj-, AjAk = 0, k, A? = Aj, A^E^ = E^A^-, 

we get the following formulas. 



< n,- >= ixr [A,S] - 1 + iyRSA,SS,Ry, (35) 
<7n,n, = ^Tr [A,SAfeS] + ^yRSA,SAfcSE,Ry, j^k (36) 
(7n,n, = ^Tr [A,-S (A,S - l2iv)] + ^yR (2SA,- - Isat) SA,-SE,Ry, (37) 

where 

CTujUk =< TT-jlT-k > - < rij >< Uk >, 

and matrix S reads 

S = (Ej;R+ l2Ar) 

6 Physical applications 

1. Suppose we have some Gaussian distribution function 

>V(q) = [det(27rM)]-^exp{-i(q- (q))M-^(q- (q))|, (38) 

where q is an iV-dimcnsional vector, and M is a symmetric positively definite N x A'^-matrix, 
which is nothing but the covariance matrix, 

M = \\Mij\\ , Mij = (qiqj) - {qi){qj), 

if id)) = I /(q)>V(q)d^q, J W(q)d^q = 1. 
In this case the characteristic function x(a) = (exp(aq)) is also Gaussian, 

X(a) = exp (^aMa + a(q)) . (39) 

On the other hand, evidently, 

oo n-i riN 

ni,...nN=0 "'l' 

Comparing these formulas with eq. (1) we arrive at the relation 

2. If A" = 2, q = (p, g), and W(q) is the Wigner function of some one-dimensional 

quantum system with quadrature components p, q, then the probability to observe n quanta 
in the state described by W(q) (in general, it is a mixed squeezed quantum state) is given by 
the relation [3,5,6,9,17] 



with the following elements of matrix R, 
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2 [(2(7pp - 1 + 2iapq) (g) + i (1 - 2(7^^ + 2iapq) (p)] 

2r - 4d - 1 



and vector y, 



Here 



vi = y2 



(i = detM > -, 
4 



T = TrM > 1 



(42) 



(43) 



(44) 



(the inequalities are the consequences of the generalized uncertainty relations [3], which 
ensure the positive definiteness of the statistical operator corresponding to function W(q)). 
The probability to have no photons Vq is given by the formula 



Po = ((i + ^T + i)-^exp 



{pfi2a,g + 1) + {q)^i2app + 1) - 4ap,(p)(g) 
1 + 2r + 4d 



(45) 



Various special cases of this photon distribution function were investigated also in [16,18]. 
Confining ourselves to the case of (p) — (q) — and taking into account eq. (4), we can 
simplify (41) as follows, 



^+2+4 



2 Ud+l- 2T" 
[4d + T+2T, 
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Ad-1 



(Ad + if - 4T2 



(46) 



For pure quantum states, when (i = 1/4, we have an oscillating photon distribution function 
[5,19,20]: 



r 2m 



2 Y (2m)! fT-1 
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(47) 



To see the "smoothing" of these oscillations for mixed quantum states we use eq. (25). Since 
in the case under study 

1-Ad 
r = — , < 0, 

2VT2 - Ad ~ 

we choose the values of T and d in the domain 

1 1 
T - - > 2d> -, 

2 ~ 2 

to ensure the inequality \r\ < 1. Then we have for n ^ 1 



2VT^-Ad + Ad-1 ' +(-1)" 2VT^-Ad-Ad+1 



Vn ~ ^ ^ ^ ' y 1 ^ . (48) 

0m(l + 2T + 4d)"+2 (r2 - Ad)^ 

For d — 1/A this expression coincides with the Stirhng asymptotics of formula (47). Eq. 
(48) shows that the oscillations of the photon distribution function disappear provided 



n 



- 4d 



> 1. 



For highly squeezed mixed states, when » 4d ^ 1, this criterion assumes the form 
2nd/T > 1. 

3. For a quantum harmonic oscillator with a time-dependent frequency the transition 
probability from the n-th to the m-th stationary level can be expressed as [1-3,7,8] 



nlml I 



with matrix R of the form 



(49) 
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and 



The explicit form of the parameters ^ and rj is determined by the solutions of the classical 
equation 

e + a;^(i)e = 0, 

moreover, the ratio may be interpreted as a reflection coefficient from some "effective 

potential barrier" a;^(t), provided a "transmitted wave" is exp(ia;t) for t — > — oo. The 
parameters Xi, X2 are expressed through the convolution of the function e{t) with the external 
time-dependent force. If this force is absent, xi = X2 = 0. Confining ourselves with this 
simplest case we can rewrite eq. (49), due to cqs. (5) and (6), in the following equivalent 
forms (remind that m and n must be either both even or both odd): 
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(51) 



This is just the case when the argument of the Jacobi and Gegenbauer polynomials 
belongs to the interval [—1, 1]. Therefore the asymptotics for iin^n ^ 1, In — m] ~ 0{1) is 
given by eq. (26): 
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If MimarctanlTyl » 1, the asymptotics of the Bessel function yields 
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(53) 



If \n — m\ ^ iinrn ~ C'(l), then eq. (28) leads to the asymptotics 
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For 1^1 ~ 0{1) this formula can be simplified as follows, 
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TTU, 



nm l^nm '1^1 
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(55) 



Two other special cases correspond to small and large values of the "reflection coefficient 
l?]/,^^. If \rj\ (this limit corresponds to an adiabatic change of the frequency a;(t), when 
1^1 — >• 1), we get from eq. (51) 



fJ'r, 



^ \n-m\ y^ 



\n—m\ 



(56) 



For large values of n and m this expression coincides with the limits of eqs. (52) and (55) 
for \ri\ 0, if one takes into account formula (22). 

An interesting phenomenon is observed when |^| ^ \ri\ 3> 1 (more precisely, when |^| ^ 
Unm)- Then for even values of n and m we get from eq. (51) 



1 (2fc-l)!!(2/-l)!! 



lei (2fc)n (2/)!! ' 
whereas the transition probability between the neighbouring odd levels is much less: 

2 



TT/2/+1 



^Ukl + 1 



w. 



21 
2k- 



(57) 



(58) 



For instance. 



So, for |,^| ^ the transition probabilities oscillate like the photon distribution function 
in eq. (47). Such a coincidence is not accidental, due to the relation [3] 



S being the normalized energy of quantum fluctuations in the final Gaussian state which has 
originated from the initial ground state with S^.^ = |. It is well known [3,21] that Gaussian 
states with large energy of fluctuations are highly squeezed, and a manifestation of squeezing 
is just the oscillating structure of various distributions related to the number of quanta. 
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